!! Copyright (C) 2009,2010,2011,2012  Marco Restelli
!!
!! This file is part of:
!!   LDGH -- Local Hybridizable Discontinuous Galerkin toolkit
!!
!! LDGH is free software: you can redistribute it and/or modify it
!! under the terms of the GNU General Public License as published by
!! the Free Software Foundation, either version 3 of the License, or
!! (at your option) any later version.
!!
!! LDGH is distributed in the hope that it will be useful, but WITHOUT
!! ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
!! or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
!! License for more details.
!!
!! You should have received a copy of the GNU General Public License
!! along with LDGH. If not, see <http://www.gnu.org/licenses/>.
!!
!! author: Marco Restelli                   <marco.restelli@gmail.com>


!>\brief
!!
!! Simplified interface to UMFPACK.
!!
!! \n
!!
!! This module provides a simplified interface to the UMFPACK solver.
!! The solution of a linear system is obtained with three function
!! calls (here \c m is the system matrix, \c x is the unknown vector
!! and \c y is the right hand side):
!! \code
!!   call umfpack_factor(m)
!!   call umfpack_solve(m,x,y)
!!   call umfpack_clean()
!! \endcode
!!
!! \note There can be an arbitrary number of calls to \c umfpack_solve
!! with the same system matrix and different right hand side. However,
!! one can not define a new system with a call to \c umfpack_factor
!! before cleaning the previous system with \c umfpack_clean.
!!
!! \note The public variable \c umfpack_print_level should be set
!! before a call to any of the module functions.
!<----------------------------------------------------------------------
module mod_umfintf

!-----------------------------------------------------------------------

 use mod_messages, only: &
   mod_messages_initialized, &
   error,   &
   warning, &
   info

 use mod_kinds, only: &
   mod_kinds_initialized, &
   wp

 use mod_sparse, only: &
   mod_sparse_initialized, &
   ! sparse types
   t_col,       &
   transpose,   &
   matmul,      &
   diag,        &
   spdiag,      &
   clear

 use mod_umfpack, only: &
   mod_umfpack_constructor, &
   mod_umfpack_destructor,  &
   mod_umfpack_initialized, &
   umfpack_control, &
   umfpack_prl,     &
   umfpack_info,    &
   umf_int,         &
   umf_dp,          &
   umf_void,        &
   umfpack_a,       &
   umfpack_at,      &
   umfpack_aat,     &
   umfpack_pt_l,    &
   umfpack_l,       &
   umfpack_lt_p,    &
   umfpack_lat_p,   &
   umfpack_lt,      &
   umfpack_lat,     &
   umfpack_u_qt,    &
   umfpack_u,       &
   umfpack_q_ut,    &
   umfpack_q_uat,   &
   umfpack_ut,      &
   umfpack_uat,     &
   umf4def,         &
   umf4pcon,        &
   umf4pinf,        &
   umf4sym,         &
   umf4num,         &
   umf4sol,         &
   umf4solr,        &
   umf4fsym,        &
   umf4fnum

!-----------------------------------------------------------------------
 
 implicit none

!-----------------------------------------------------------------------

! Module interface

 public :: &
   mod_umfintf_constructor, &
   mod_umfintf_destructor,  &
   mod_umfintf_initialized, &
   umfpack_print_level, &
   umfpack_factor, &
   umfpack_solve,  &
   umfpack_clean

 private

!-----------------------------------------------------------------------

 ! UMFPACK related variables
 integer, target :: umfpack_print_level
 ! the two following variables are used to keep track of status of
 ! the two symbolic and numeric umfpack objects
 logical :: symbolic_set, numeric_set
 integer(umf_void) :: symbolic, numeric
 real(umf_dp) :: umf_contr(umfpack_control), umf_info(umfpack_info)

 logical, protected :: &
   mod_umfintf_initialized = .false.
 character(len=*), parameter :: &
   this_mod_name = 'mod_umfintf'

!-----------------------------------------------------------------------

contains

!-----------------------------------------------------------------------

 subroutine mod_umfintf_constructor()
  character(len=*), parameter :: &
    this_sub_name = 'constructor'

   !Consistency checks ---------------------------
   if( (mod_messages_initialized.eqv..false.) .or. &
          (mod_kinds_initialized.eqv..false.) .or. &
         (mod_sparse_initialized.eqv..false.) ) then
     call error(this_sub_name,this_mod_name, &
                'Not all the required modules are initialized.')
   endif
   if(mod_umfintf_initialized.eqv..true.) then
     call warning(this_sub_name,this_mod_name, &
                  'Module is already initialized.')
   endif
   !----------------------------------------------

   call mod_umfpack_constructor()
   symbolic_set = .false.
   numeric_set  = .false.

   mod_umfintf_initialized = .true.
 end subroutine mod_umfintf_constructor

!-----------------------------------------------------------------------
 
 subroutine mod_umfintf_destructor()
  character(len=*), parameter :: &
    this_sub_name = 'destructor'
   
   !Consistency checks ---------------------------
   if(mod_umfintf_initialized.eqv..false.) then
     call error(this_sub_name,this_mod_name, &
                'This module is not initialized.')
   endif
   !----------------------------------------------

   call umfpack_clean()
   call mod_umfpack_destructor()

   mod_umfintf_initialized = .false.
 end subroutine mod_umfintf_destructor

!-----------------------------------------------------------------------

 subroutine umfpack_factor(m)
  type(t_col), intent(in) :: m
 
  character(len=*), parameter :: &
    this_sub_name = 'umfpack_factor'
 
   ! Inizialize umf_contr with the default UMFPACK parameters
   call umf4def(umf_contr)
   ! set the verbose mode
   umf_contr(umfpack_prl) = real( umfpack_print_level , umf_dp )
   call umf4pcon(umf_contr)          ! output
   ! pre-order and symbolic analysis
   if(symbolic_set) call error(this_sub_name,this_mod_name, &
   'Trying to factor a system when symbolic is already set.')
   call umf4sym( int(m%m,umf_int) , int(m%n,umf_int) ,  &
                int(m%ap,umf_int) , int(m%ai,umf_int) , &
                real(m%ax,umf_dp) , symbolic,           &
                        umf_contr , umf_info ); symbolic_set = .true.
   call umf4pinf(umf_contr,umf_info) ! output
   ! numeric factorization
   if( numeric_set) call error(this_sub_name,this_mod_name, &
   'Trying to factor a system when numeric is already set.')
   call umf4num( int(m%ap,umf_int) , int(m%ai,umf_int) , &
                 real(m%ax,umf_dp) , symbolic,numeric ,  &
                         umf_contr , umf_info); numeric_set = .true.
   call umf4pinf(umf_contr,umf_info) ! output
   ! free the symbolic object
   call umf4fsym(symbolic); symbolic_set = .false.

 end subroutine umfpack_factor
 
!-----------------------------------------------------------------------
 
 subroutine umfpack_solve(m,x,y,transposed_mat)
  type(t_col), intent(in) :: m
  real(wp),    intent(in) :: y(:)
  real(wp),   intent(out) :: x(:)
  logical, optional, intent(in) :: transposed_mat
 
  integer(umf_int) :: system_type
  real(umf_dp), allocatable :: x_dp(:)
  character(len=*), parameter :: &
    this_sub_name = 'umfpack_solve'
 
   if(.not.numeric_set) call error(this_sub_name,this_mod_name, &
     'Can not solve a system if "numeric" is not set.')

   system_type = umfpack_a
   if(present(transposed_mat)) then
     if(transposed_mat) system_type = umfpack_aat
   endif

   ! solve
   allocate(x_dp(size(x)))
   call umf4solr(system_type,int(m%ap,umf_int),int(m%ai,umf_int), &
     real(m%ax,umf_dp),x_dp,real(y,umf_dp),numeric,umf_contr,umf_info)
   call umf4pinf(umf_contr,umf_info) ! output
   x = real(x_dp,wp)
   deallocate(x_dp)

 end subroutine umfpack_solve
 
!-----------------------------------------------------------------------
 
 subroutine umfpack_clean()
  character(len=*), parameter :: &
    this_sub_name = 'umfpack_clean'
 
   if(symbolic_set) then
     call umf4fsym(symbolic); symbolic_set = .false.
   endif
   if(numeric_set ) then
     call umf4fnum(numeric ); numeric_set  = .false.
   endif

 end subroutine umfpack_clean
 
!-----------------------------------------------------------------------

end module mod_umfintf

